#ifndef __EVAL_THRESH_H__
#define __EVAL_THRESH_H__

#include "Image.H"

/**
 * Computes Accurancy, Sensitivity, Specificity and Overlap.
 * Is templated
 *
 * @param reference   Input image
 * @param result      Input image
 * @param accuracy    Output statistical-value
 * @param sensitivity Output statistical-value
 * @param specificity Output statistical-value
 * @param overlap	  Output statistical-value
 * @param thres       Input Threshold value
 */
template <typename T, typename U>
void evaluate(const pip::Image<T> &reference, 
			  const pip::Image<U> &result, 
              float &accuracy, 
              float &sensitivity, 
              float &specificity, 
              float &overlap,
			  float const &thres=0.5)
{
  int TP, TN, FP, FN;					//init values
  TP = TN = FP = FN = 0;
  for (int z = 0; z < result.dim(2); ++z) {
	for (int y = 0; y < result.dim(1); ++y) {
	  for (int x = 0; x < result.dim(0); ++x) { 
		if(float(result(x,y,z)) > thres)		 //for all pixel values check if they are greater than threshold 
		 {
			  if(reference(x,y,z)==1) TP++;	    // if they are and the reference is positive - IT is a true  positive
			  else FP++;							                           //is negative - IT is a false positive
		 }else{
			  if(reference(x,y,z)==1) FN++;     // if not and the reference is positive - IT is a false negative
			  else TN++;												  //is negative - IT is a true  negative
		 }
	   }
     }	
  }

 accuracy    = float((TP + TN ) /float(TP+TN+FP+FN));  // equations according to the wikipage
 sensitivity = float( TP        /float(TP+FN));
 specificity = float( TN        /float(TN+FP));
 overlap     = float( TP        /float(TP+FP+FN));

};

/**
 * Performs optimal thresholding. Computes a segmentation of the input image based
 * on an optimal thresholding approach, as explained in Sonka, algorithm 6.2.
 *
 * @param input  Input  image
 * @param output Output image
 * @param threshold Threshold value
 * @param mask   Input  image-mask
 */
void optimal_thresholding(const pip::Image<short> &input,
						  pip::Image<unsigned char> &output,
                          float &threshold,
						  const pip::Image<unsigned char> &mask = pip::Image<unsigned char>())
{
	//if input and output are of equal size AND mask is empty OR of equal size start processing
	if((input.dimension()==output.dimension()) && ( mask.dim(0)==0 || (input.dimension()==mask.dimension()) ))
	{
		float first_thresh = 0;
		if(mask.dim(0)==0)   //if no mask is passed use corner-pixel for first threshold
		{
			//get sum of all pixels
			float sum_all(0), sum_corner(0);
			int	  num_all(0), num_corner(4);
			for (int z = 0; z < input.dim(2); ++z) {
				for (int y = 0; y < input.dim(1); ++y) {
					for (int x = 0; x < input.dim(0); ++x) {
						sum_all += float(input(x,y,z));  // count object pixel and sum their values 
						++num_all;
					}
				}
			}
			//get sum of corner pixel
			sum_corner = float(
								 input(0,0,0)                      //initial threshold determined by the corner-pixels
								+input(input.dim(0)-1,0,0)				//for 3d - only in the 0-slice 
								+input(0,input.dim(1)-1,0)
								+input(input.dim(0)-1,input.dim(1)-1,0));
			
			if(input.dim(2)>0) //if 3d one needs more corner pixel
			{
				sum_corner   += input(0,0                          ,input.dim(2)-1) //next pixels for initial threshold
							   +input(input.dim(0)-1,0             ,input.dim(2)-1) //for 3d - only in the dim(2)-slice 
							   +input(0,input.dim(1)-1             ,input.dim(2)-1)
							   +input(input.dim(0)-1,input.dim(1)-1,input.dim(2)-1);
				num_corner   += 4;
			}
			//calculate first threshold and check for DEVISION by ZEROS 
			if(!num_all)  // if only 4 pixels in the image (or 8 in 3d)
			{
				first_thresh = (sum_corner/num_corner)/2; 
			}
			else
			{
				first_thresh = (sum_corner/num_corner + ( (sum_all-sum_corner)/(num_all-num_corner)) )/2; 
			}

		}
		else    // if mask is passed - we already checked if it is correct size
		{		// use min and max values outside mask for first threshold
			float t_min(9999);
			float t_max(-9999);
			for (int z = 0; z < input.dim(2); ++z) {
				for (int y = 0; y < input.dim(1); ++y) {
					for (int x = 0; x < input.dim(0); ++x) {
						if(int(mask(x,y,z)) > 0)				//if not in mask
						{
							if(t_min > float(input(x,y,z)))			//find and remember min and max
							{
								t_min = float(input(x,y,z));
							}
							else if(t_max < float(input(x,y,z)))
							{
								t_max = float(input(x,y,z));
							}
						}
					}
				}
			}
			first_thresh = (t_min+t_max)/2;
		}

		std::cout << "Initial threshold: " << first_thresh << std::endl;  // print the first threshold used

		float last_thresh(-9999);
		float sum_pos(0);
		float sum_neg(0);
		float num_pos(0);
		float num_neg(0);
		float mean_obj(-9999);
		float mean_bgr(-9999);

		do{				// iterative part !!! 
			sum_pos = sum_neg = num_pos = num_neg = 0;
			for (int z = 0; z < input.dim(2); ++z) {
				for (int y = 0; y < input.dim(1); ++y) {
					for (int x = 0; x < input.dim(0); ++x) {
						if(mask.dim(0) == 0 || mask(x,y,z)){
							if(float(input(x,y,z)) > first_thresh)  //use threshold to...
							{
								sum_pos += float(input(x,y,z));  // count object pixel and sum their values 
								++num_pos;
							}
							else
							{
								sum_neg += float(input(x,y,z));  // count background pixel and sum their values 
								++num_neg;
							}
						}
					}
				}
			}
			if(num_pos)			//check for devision by zero
			{
				mean_obj = sum_pos/num_pos;   // get mean values
			}
			else
			{
				mean_obj = 0;
			}
			
			if(num_neg)			//check for devision by zero
			{
				mean_bgr = sum_neg/num_neg;   // get mean values
			}
			else
			{
				mean_bgr = 0;
			}
			last_thresh  = first_thresh;      // remember last threshold
			first_thresh = (mean_obj+mean_bgr)/2;  // and get new threshold
			std::cout << "Threshold: " << first_thresh << std::endl;
		
		}while(first_thresh!=last_thresh);   // test if any progress has been made in the mean time
		// alternatively use  while( (first_thresh-last_thresh) < epsilon)
		// to avoid problems because of numeric precision...
		// it wasn't necassary for the integer valued images used in this assignment


		for (int z = 0; z < input.dim(2); ++z) {
			for (int y = 0; y < input.dim(1); ++y) {
				for (int x = 0; x < input.dim(0); ++x) {
					if(float(input(x,y,z)) > first_thresh)		// create the binary[0,1] output image
					{
						output(x,y,z)=1;	
					}else{
						output(x,y,z)=0;
					}
				}
			}
		}
	}
	else
	{
		pip::error("", "The two Images are not compatible. Only use equal sized images!\n");
 		return;
	}	
};



#endif 